function y=two_regimes(xxd)
clc;

rho=0.05;
mu_L=0.5;
mu_H=0.5;
eta_1=0.3;
eta_2=0.3;
phi_H1=2*2;
phi_L1=-1*2;
phi_H2=phi_H1+xxd;
phi_L2=phi_L1-xxd;

% No amplitude change
BB=[(rho+mu_L)/phi_H1 -mu_L/(rho+mu_L)*(rho+mu_L)/phi_H1;-mu_H/(rho+mu_H)*(rho+mu_H)/phi_L1 (rho+mu_H)/phi_L1];

dd=eig(BB);
om1=dd(2);
om2=dd(1);
Z=om1^2/om2^2 *(1-om2*phi_H1/(rho+mu_L))/(1-om1*phi_H1/(rho+mu_L));
format long
tt=1/(om2-om1)*log(Z)


BB2=[(rho+mu_L)/phi_H2 -mu_L/(rho+mu_L)*(rho+mu_L)/phi_H2;-mu_H/(rho+mu_H)*(rho+mu_H)/phi_L2 (rho+mu_H)/phi_L2];

dd2=eig(BB2);
Om1=dd2(2);
Om2=dd2(1);
Z2=Om1^2/Om2^2 *(1-Om2*phi_H2/(rho+mu_L))/(1-Om1*phi_H2/(rho+mu_L));
format long
tt2=max(1/(Om2-Om1)*log(Z2),0)






% With amplitude change
% Code
A1=[rho+mu_L+eta_2,-eta_2,-mu_L,0];
A2=[-eta_1,rho+mu_L+eta_1,0,-mu_L];
A3=[-mu_H,0,rho+mu_H+eta_2,-eta_2];
A4=[0,-mu_H,-eta_1,rho+mu_H+eta_1];

A=[A1;A2;A3;A4];
M=diag([1/phi_H1,1/phi_H2,1/phi_L1,1/phi_L2]);
B=M*A;
[Q,Omega]=eig(B);

om_1=Omega(1,1);
om_2=Omega(2,2);
om_3=Omega(3,3);
om_4=Omega(4,4);


Z1=@(om_st) [om_1*exp(om_1*om_st(1)), om_2*exp(om_2*om_st(1)), om_3*exp(om_3*om_st(1)), om_4*exp(om_4*om_st(1));
             om_1*exp(om_1*om_st(2)), om_2*exp(om_2*om_st(2)), om_3*exp(om_3*om_st(2)), om_4*exp(om_4*om_st(2));
            1,1,1,1;
            1,1,1,1];

Z=@(om_st) Z1(om_st).*Q;

C=@(om_st) inv(Z(om_st))*[1;1;0;0];

BBB=@(om_st) [om_1^2*exp(om_1*om_st(1)), om_2^2*exp(om_2*om_st(1)), om_3^2*exp(om_3*om_st(1)), om_4^2*exp(om_4*om_st(1));
              om_1^2*exp(om_1*om_st(2)), om_2^2*exp(om_2*om_st(2)), om_3^2*exp(om_3*om_st(2)), om_4^2*exp(om_4*om_st(2))];

d=@(om_st) (C(om_st))'*(BBB(om_st).*[Q(1,:);Q(2,:)])';
                                    

f=@(om_st) sqrt(d(om_st)*d(om_st)');

options = optimset('PlotFcns','optimplotfval','TolX',1e-16, 'TolFun',1e-16);

[X,FVAL]=fminsearch(f,[tt tt2],options)

y=[tt,tt2,X(1),X(2)]
return
